PrecipitationRead Subroutine

SUBROUTINE PrecipitationRead &
!
( time, dem )

IMPLICIT NONE

!Arguments with intent(in):
TYPE (DateTime), INTENT(IN) :: time !!current time
TYPE (grid_real), INTENT(IN) :: dem

!local declarations:
TYPE (DateTime) :: time_toread !! time to start reading from
CHARACTER (LEN = 300) :: filename
CHARACTER (LEN = 300) :: varname
INTEGER (KIND = short) :: i, j

!-------------------------end of declarations----------------------------------

IF ( .NOT. (time < timeNew) ) THEN
   
   !time_toread = time + - (dtPrecipitation - raingages % timeIncrement)
   time_toread = time 
   timeString = time_toread
   CALL Catch ('info', 'Precipitation',   &
		                 'read new precipitation data: ', &
                     argument = timeString)
   
   !update lapse rate when required
    IF (elevationDrift == 1) THEN
      IF (time == lapse_map % next_time) THEN !update lapse_map
          varname = lapse_map %var_name
          filename = lapse_map % file_name
          !destroy current grid
          CALL GridDestroy (lapse_map)
          !read new grid in netcdf file
          CALL NewGrid (layer = lapse_map, fileName = TRIM(filename), &
                         fileFormat = NET_CDF, &
                         variable = TRIM(varname), &
                         time = time_toread)
      END IF
    END IF

   SELECT CASE (interpolationMethod)
    CASE (0) !input grid

		CALL ReadField (fileGrid,  time_toread, &
        dtPrecipitation, dtGrid, &
        'C', precipitationRate, &
            varName = precipitationRate % var_name)

      
    CASE DEFAULT !requires interpolation
      !read new precipitation data
      CALL ReadData (network = raingages, fileunit = fileunit, &
                      time = time_toread, aggr_time = dtPrecipitation, &
                      aggr_type = 'sum', tresh = valid_prcn)
      
      IF (elevationDrift == 1) THEN

        !shift precipitation observation to reference elevation by applying lapse rate
        CALL ShiftMeteoWithLapse (raingages, lapse_map, refElevation, &
                                  stationsRefElev, dtPrecipitation)
        
        !interpolate 
        IF (interpolationMethod_assignment == 1) THEN !only one method is applied
                
                CALL Interpolate (network = stationsRefElev, &
                            method = interpolationMethod, &
                            near = neighbors, &
                            idw_power = idw_power, & 
                            anisotropy = krige_anisotropy, &
                            varmodel = krige_varmodel, &
                            lags = krige_lags, &
                            maxlag = krige_maxlag, &
                            grid = precipitationRate, &
                            devst = grid_devst)                  
                
        ELSE 
            !loop trough interpolation methods
            DO i = 1, 3
                IF (interpolationMethod_vector(i) > 0) THEN
                      
                    CALL Interpolate (network = stationsRefElev, &
                            method = interpolationMethod_vector(i), &
                            near = neighbors, &
                            idw_power = idw_power, & 
                            anisotropy = krige_anisotropy, &
                            varmodel = krige_varmodel, &
                            lags = krige_lags, &
                            maxlag = krige_maxlag, &
                            grid = interpolatedMap (i), &
                            devst = grid_devst) 
                      
                END IF
            END DO

            !compose the final interpolated field
            DO i = 1, interpolationMethod_map % idim
                DO j = 1, interpolationMethod_map % jdim
                    IF (interpolationMethod_map % mat(i,j) /= &
                        interpolationMethod_map % nodata ) THEN
                        precipitationRate % mat (i,j) = &
                        interpolatedMap (interpolationMethod_map % mat(i,j)) % mat(i,j)
                    END IF
                END DO
            END DO
        END IF

        !Shift back interpolated field to terrain surface
        CALL ShiftBackWithLapse (precipitationRate, dem, &
                                lapse_map, refElevation, &
                                dtPrecipitation)

      ELSE
         
        IF (interpolationMethod_assignment == 1) THEN !only one method is applied
                
                CALL Interpolate (network = raingages, &
                            method = interpolationMethod, &
                            near = neighbors, &
                            idw_power = idw_power, & 
                            anisotropy = krige_anisotropy, &
                            varmodel = krige_varmodel, &
                            lags = krige_lags, &
                            maxlag = krige_maxlag, &
                            grid = precipitationRate, &
                            devst = grid_devst) 
                
        ELSE
            !loop trough interpolation methods
            DO i = 1, 3
                IF (interpolationMethod_vector(i) > 0) THEN
                      
                    CALL Interpolate (network = raingages, &
                            method = interpolationMethod_vector(i), &
                            near = neighbors, &
                            idw_power = idw_power, & 
                            anisotropy = krige_anisotropy, &
                            varmodel = krige_varmodel, &
                            lags = krige_lags, &
                            maxlag = krige_maxlag, &
                            grid =  interpolatedMap (i), &
                            devst = grid_devst) 
                    
                END IF
            END DO

            !compose the final interpolated field
            DO i = 1, interpolationMethod_map % idim
                DO j = 1, interpolationMethod_map % jdim
                    IF (interpolationMethod_map % mat(i,j) /= &
                        interpolationMethod_map % nodata ) THEN
                        precipitationRate % mat (i,j) = &
                        interpolatedMap (interpolationMethod_map % mat(i,j)) % mat(i,j)
                    END IF
                END DO
            END DO

        END IF  !1 or more interpolation methods       
      END IF  !elevationDrift
		
	END SELECT

  !apply scale factor and offset, if defined
	IF (offset_value /= MISSING_DEF_REAL) THEN
	   CALL Offset (precipitationRate, offset_value)
	END IF
	
	
	IF (scale_factor /= MISSING_DEF_REAL) THEN
	   CALL Scale (precipitationRate, scale_factor)
    END IF

    
    !filter data < 0.
    DO i = 1, precipitationRate % idim
      DO j = 1, precipitationRate % jdim
           IF (precipitationRate % mat(i,j) /= precipitationRate % nodata) THEN
              IF ( precipitationRate % mat(i,j) < 0.) precipitationRate % mat(i,j) = 0.
           END IF
      END DO
    END DO

    

  !grid exporting
  IF(export > 0 ) THEN
	  IF( time == timeNewExport .AND. &
        time >= export_start .AND. &
        time <= export_stop ) THEN
        timeString = time
        timeString = timeString(1:19)
        timeString(14:14) = '-'
		    timeString(17:17) = '-'
        
        grid_devst % reference_time = precipitationRate % reference_time
        IF (needConversion) THEN
           CALL GridConvert (precipitationRate, exportedGrid)
           CALL GridConvert (grid_devst, exportedGridVar)
        ELSE
           exportedGrid = precipitationRate
           exportedGridVar = grid_devst
        END IF 

        SELECT CASE (export_format)
        CASE (1) !esri-ascii
              export_file = TRIM(export_path) //  TRIM(timeString) // &
                           '_precipitation.asc'
              CALL Catch ('info', 'Precipitation',   &
		                       'exporting grid to file: ', &
                           argument = TRIM(export_file))
              CALL ExportGrid (exportedGrid, export_file, ESRI_ASCII)
              
              IF (krige_var == 1) THEN !export kriging variance
                    export_file_var = TRIM(export_path) //  TRIM(timeString) // &
                           '_precipitation_variance.asc'
                    CALL Catch ('info', 'Precipitation',   &
		                       'exporting variance grid to file: ', &
                              argument = TRIM(export_file_var))
                    CALL ExportGrid (exportedGridVar, export_file_var, ESRI_ASCII)
              END IF
              
        CASE (2) !esri-binary
              export_file = TRIM(export_path) //  TRIM(timeString) // &
                           '_precipitation'
              CALL Catch ('info', 'Precipitation',   &
		                       'exporting grid to file: ', &
                           argument = TRIM(export_file))
              CALL ExportGrid (exportedGrid, export_file, ESRI_BINARY)
              
              IF (krige_var == 1) THEN !export kriging variance
                   export_file_var = TRIM(export_path) //  TRIM(timeString) // &
                               '_precipitation_variance'
                   CALL Catch ('info', 'Precipitation',   &
		                       'exporting variance grid to file: ', &
                              argument = TRIM(export_file_var))
                   CALL ExportGrid (exportedGridVar, export_file_var, ESRI_BINARY)
              END IF
              
        CASE (3) !net_cdf
              CALL SetCurrentTime (time, exportedGrid)
              CALL Catch ('info', 'Precipitation',   &
                           'exporting grid to file: ', &
                           argument = TRIM(export_file))
              CALL ExportGrid (exportedGrid, export_file, 'precipitation_amount', 'append')
              
              IF (krige_var == 1) THEN !export kriging variance
                  CALL SetCurrentTime (time, exportedGridVar)
                  CALL Catch ('info', 'Precipitation',   &
		                       'exporting grid to file: ', &
                              argument = TRIM(export_file_var))
                 CALL ExportGrid (exportedGridVar, export_file_var, 'precipitation_amount_variance', 'append')
              END IF
        END SELECT
        timeNewExport = timeNewExport + export_dt
    END IF
  ENDIF


   
  !conversion mm => m/s 
  DO i = 1, precipitationRate % idim
    DO j = 1, precipitationRate % jdim
       IF (precipitationRate % mat(i,j) /= precipitationRate % nodata) THEN
           precipitationRate % mat(i,j) = precipitationRate % mat(i,j) / &
           dtPrecipitation / 1000.
       END IF
    END DO
  END DO


  !time forward
  timeNew = timeNew + dtPrecipitation
END IF

RETURN
END SUBROUTINE PrecipitationRead